# source functions

pred_prob_no_het = function(beta0,beta1,epsilon,x){
  ystarmean = beta0 + beta1*x + mean(epsilon)
  ystar25 = beta0 + beta1*x + quantile(epsilon,probs=0.025)
  ystar975 = beta0 + beta1*x + quantile(epsilon,probs=0.975)
  pr_mean = exp(ystarmean)/(1+exp(ystarmean))
  pr_25 = exp(ystar25)/(1+exp(ystar25))
  pr_975 = exp(ystar975)/(1+exp(ystar975))
  plot(NA,xlim = c(min(x),max(x)),ylim=c(0,1),
       xlab = "X value",
       ylab = "Predicted probability")
  lines(x,pr_mean)
  lines(x,pr_25,lty="dashed")
  lines(x,pr_975,lty="dashed")
}

dald = function(x,p){
  out = ifelse(x<0,p*(1-p)*exp(-x*(p-1)),p*(1-p)*exp(-x*p) )
  return(out)
}

pald = function(x,p){
  out = ifelse(x<0,p*exp(x*(1-p)), 1-(1-p)*exp(-x*(p)))
  return(out)
}

inverse = function (f,p, lower = -10000, upper = 10000) {
  function (y,p) uniroot((function (x) f(x,p) - y), lower = lower, upper = upper)[1]
}

qald = inverse(pald)

rald = function(n,p){
  tmp = runif(n)
  out =  NA
  for (i in 1:n){
    out[i] = qald(tmp[i],p)
  }
  return(unlist(out))
}
